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A simple model for simulating the current dynamics and the magnetic properties of superconduc- 
tors is presented. Short simulation times are achieved by solving the differential form of Maxwell's 
equations inside the sample, whereas integration is only required at the surface to meet the exact 
boundary conditions. The procedure reveals the time and position dependence of the current den- 
sity and the magnetic induction (B) making it very convenient to apply field dependent material 
parameters for the simulation of magnetization loops, relaxation measurements, etc. Two examples, 
which are important for standard magnetometry experiments, are discussed. Firstly, we prove that 
evaluating the critical current density ( J c ) from experiment by applying Bean's model reveals almost 
the exact J C (B) behavior, if the evaluation is corrected by a simple numerical expression. Secondly, 
we show that the superconducting volume fraction of a sample can be directly determined from 
magnetization loops by carefully comparing experiment and simulation in the field range, where the 
current loops are differently oriented within the sample. 
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I. INTRODUCTION 

Simulating the magnetic behavior of type II super- 
conductors requires solving Maxwell's equations under 
specific material equations. Many different models have 
been proposed for solving this task including energy min- 
imization, variational, or more direct approaches (e.g. 
Refs. EyiimHSIl). When addressing not only the 
static, but also the dynamic behavior, it is often con- 
venient to calculate the time evolution of field relevant 
quantities (e.g. of the magnetic vector potential or of the 
induction) directly. 

Still, different approaches can be chosen. The most 
direct way is to solve Maxwell's equations by integra- 
tion (i.e. with the help of the appropriate Green's func- 
tion) which yields, e.g., the current density or the vector 
potential^ This method naturally leads to results fulfill- 
ing the electromagnetic boundary conditions at the sam- 
ple surface. However, integrating over the whole sample 
volume requires a lot of calculation time and is usually 
too slow for three-dimensional (3D) problems. For a grid 
of iV 3 elements, into which the sample is divided for the 
numerics, at least TV 6 operations have to be performed 
per time step. 

The opposite approach deals with the differential form 
of Maxwell's equations, where fields and currents are lo- 
cally related by difference equations^ This makes the 
calculations very fast inside the sample, since only about 
N 3 operations per time step are required. Care has to 
be taken, however, with the boundary conditions at the 
surface, which requires solving the equations also outside 
the sample up to distances, where the influence of the 
sample can be neglected. 

As a compromise, we present a model, where the dif- 
ferential equations are applied to the sample interior, 
whereas the values at the sample boundaries are obtained 
by integrating over the sample volume yielding the cor- 



rect boundary conditions. This requires N 3 + 6N 4 op- 
erations. It was already demonstrated^ that the model 
worked sufficiently rapidly for 3D calculations on a stan- 
dard PC. 

In this paper, we provide details of the method and 
discuss examples relevant for every-day magnetometry 
measurements. In section [TT] the simulation procedure 
and some technical aspects of the implementation are de- 
scribed. We start with a general model, but concentrate 
in the following on rectangular samples, for which several 
tests were carried out. Section Hill is devoted to applica- 
tions of the method to magnetometry (e.g. SQUID or 
VSM) measurements. We will verify the well known ex- 
perimental evaluation of the critical current density from 
magnetization loops by applying Bean's model and intro- 
duce a simple way of determining the superconducting 
volume fraction of a sample. 



II. SIMULATION MODEL 

Our aim is to solve Maxwell's equations relevant 
for the magnetic properties of a superconductor. We 
assume superconductors without (non-superconducting) 
para- or diamagnetic properties and initially ignore the 
reversible magnetization of the superconducting state. 
Thus, Maxwell's equations (in SI units) for vacuum may 
be applied. 



V-B = 
V x B — no J 
d t B = -V x E 



(1) 
(2) 
(3) 



The displacement current (Eq. \2§ is not considered. B 
denotes the magnetic induction, (iq — 47r x ICU 7 TmA -1 
the vacuum permeability, J the current density and E 
the electric field. Generally, the vector quantities depend 
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on time (t) and position (r). The electric field is usually 
explicitly provided, e.g. as a function of current den- 
sity, induction, and material specific parameters such as 
the critical current density to specify the superconduct- 
ing properties of the material. It is directed parallel to 
B xv, where v denotes the velocity of the vortices caused 
by the Lorentz force J x B. In the following, we will 
address rather simple configurations, relevant for most 
magnetometry experiments, where the currents flow per- 
pendicular to B and thus E and J are parallel. The more 
complicated situation, where J has also a (usually small) 
component parallel to B, which can occur in the case of 
arbitrary field directions or irregular sample geometries, 
etc., is not considered. Recently, this problem was solved 
for simple configurations in the (static) Bean limi t 11 ' 12 , 
but it is currently not establish how to obtain more gen- 
eral results. For the simpler examples, discussed here, 
different kinds of material equations are available in the 
literature, which describe different materials and experi- 
mental situations (e.g. Refs. Il3ul4ll5l ). Most prominent 
is the power law 



(4) 



which reproduces the flux creep behavior of type II su- 
perconductors quite well. Here, E c is the electric field cri- 
terion, J c the critical current density, and n the so-called 
n- value, which is a two-digit number in most cases. Very 
large n values (e.g. n — 00) reproduce the Bean models 
(E = for J < J c and E = oo for J > J c ). Both J c and 
n may depend on B and other parameters. 



A. Procedure 

i. Our starting point is a convenient initial state of the 
magnetic induction B(f,t) and the corresponding 
current density J(r,t) (satisfying V x B = ^J). 
Usually, we start with B = and J = in the 
whole sample or with the results of a previous sim- 
ulation. 

ii. Next, we calculate the electric field - E(r,t) - by 
applying the material equation (e.g. Eq. 0| and 

iii. its curl, which corresponds to the time derivative 
of the induction (d t -B, cf. Eq. [3]) and thus leads to 

iv. the induction of the next time step 

B(f,t + At) = B{f,t) + d t B(r,t)At (5) 

where At denotes a small time increment. Equa- 
tions (O and (O are only applied to the interior of 
the sample, whereas 
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FIG. 1: (Color online) Flow diagram of the simulation se- 
quence when only the irreversible magnetization is taken into 
account (see the text for details). 



v. B(r,t) at the sample surface is obtained by inte- 
grating over the sample volume to fulfill the bound- 
ary conditions 

B(r, t + At) = Mo[ffa(r, t + At) + H s (r, t)] (6) 
with 



1 

4tt 



(7) 



H a is the applied field. Equation ([6]) is the only 
place where this quantity appears (i.e. H a influ- 
ences our calculations only via the surface). 

vi. Knowing B in the whole sample allows evaluating 
J{r,t + At) fromEq. ©. Note that B(r,t + At) © 
is calculated from H s (r , t) , which actually refers to 
the current distribution of the previous time step. 
This calculation could be repeated with the new 
currents J(r, t + At), but this was found to be un- 
necessary due to the small value of At. 

vii. We end up with a new set of B and J, close our 
simulation loop and proceed with calculating the 
next E (step ii). 

Figure Q] outlines the major parts of the sequence. Af- 
ter the last step (vii) specific quantities may be evaluated 
(which is usually not done after each simulation loop, but 
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after predefined field or time intervals), e.g. the magnetic 
moment 

m = ^ J tfrfx J (8) 

Thanks to the use of the difference equations, the simu- 
lations can be carried out quite rapidly. Time-consuming 
integration over the whole sample volume is only neces- 
sary for the induction at the sample surface to provide 
the correct boundary values. 

We wish to point out another very attractive feature 
of this method, namely the fact, that the induction B 
is directly assessed in the simulations. This makes it 
very convenient to apply any field dependent property 
- in particular a field dependent critical current density 
J C (B) - without the need of any additional calculations 
of B (as would be necessary for many other approaches). 

Different specifications of E and H a allow dealing with 
different materials and experiments. For instance, the 
simple power law of Eq. (fj| leads to the typical flux creep 
behavior of superconductors. A transport current can be 
considered by adding external electric fields. 

So far, only the irreversible properties have been ad- 
dressed in the simulation procedure. A simple way of 
taking also the reversible magnetization of a supercon- 
ductor into account is presented in Ref. y, which can be 
easily adapted in our method by adding two further steps 
to the above procedure. 

viii. After finishing the last step (vii) of the original pro- 
cedure, we know B and J in the entire sample. 
B is the overall local magnetic induction and thus 
includes also the contribution from the reversible 
magnetization M T , i.e. 

B = H {H + M t ). (9) 

The reversible behavior may be assessed from 
Ginzburg Landau theory or from simpler approxi- 
mations (like those of Ref. [13), which provide M r as 
a function of B (or H) and two parameters such as 
B C 2 - the upper critical field - and k - the Ginzburg 
Landau parameter. Note that M r vanishes at the 
sample surface. 

ix. Equation is still valid and the resulting J needed 
for calculating the boundary values (Eq. ([7])). Note 
that J includes also the reversible part of the cur- 
rents and thus may be particularly large at the sam- 
ple surface. However, J does not enter the Lorentz 
force density (/l) that drives the flux lines, bul4 
/i=/hXB, with 

V x H = J H , (fO) 

where H is obtained from Eq. (|9|). Accordingly, 
Jh (but not J) generates the electric field and we 
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FIG. 2: (Color online) Simulation of the magnetization (solid 
line, M = m/volume) of a superconductor (axfexc=lx 
I x I mm 3 ), where both the reversible and the irreversible 
contributions are significant. J C (B) is taken from Eq. (|24p 
with J c0 = 2 x 10 s Am -2 , B = 0.02 T, and a = 0.5. The 
dashed line illustrates the pure reversible fraction - M r (_ff a ) 
(obtained from the approximate equations of Ref. [l7h. where 
B c2 = 1 T and k = 4.8 - thus B cl ~ 0.05 T. 

have to replace J by Jjj in the expression for E 
(i.e. in the material equation). If the reversible 
magnetization is ignored (as, e.g., in the flow chart 
of Fig. UJI, M r = and thus J — Jh. 

Step (ix) finishes the procedure, in case the reversible 
magnetization is not ignored, and we can proceed with 
step (ii). Figure [2] presents an example of a magnetiza- 
tion loop, where both the reversible and irreversible con- 
tributions are significant. For describing flux line motion 
we employ 

as suggested in Ref. |3j, which correctly leads to flux flow 
(B k BJh) at Jh » J c and flux creep at Jh <C J c - 
Here, A c is an appropriate constant, equivalent to E c of 
Eq. (JH, but given in units of VrnT" 1 A -1 . 

In the following, we mainly concentrate on samples 
with high critical current density, where the reversible 
properties may be ignored and therefore only steps (i - 
vii) will be executed. 

B. Technical aspects 

In this section, we address some specific points of our 
implementation. In most cases, we deal with rectangular 
samples and, therefore, subdivide the sample into small 
rectangular elements, within which the current density is 
kept constant. For simplicity, we work with an equidis- 
tant grid, but note that a better efficiency is expected 
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FIG. 3: (Color online) Schematic view of a single element, (in 
which J is constant) and the positions, where B, E, and J are 
calculated. The dotted square indicates the two dimensional 
case (see the text for details). 



ing the initial stage, and can be stored for later use in 

Eqs. nm, 

Finally, we address the problem of choosing an ade- 
quate time increment At, which should be as large as 
possible to make the simulation time short, but small 
enough to have no influence on the results. In our imple- 
mentation, it turns out that choosing At too large makes 
the calculations unstable with |J| increasing rapidly to- 
wards infinity. This situation can be easily monitored 
and mitigated by resuming the simulation from a previ- 
ously saved configuration with a smaller At. It turns out 
that even smaller values of At do not modify the results. 



from a grid with smaller elements near the surface, since 
most quantities typically exhibit larger variations there. 
Each grid element consists of eight vertices and a cen- 
tral point (Fig. [3]). The sequence (Fig. [T]) dictates the 
position, on which each quantity is evaluated, in a nat- 
ural way (Fig. [3]). The induction B is defined at the 
vertices, which also include the sample surface. Accord- 
ingly, its derivative (or actually its curl) is calculated at 
the midpoint position using B values of the next eight 
adjacent vertices. The electric field is in most cases a di- 
rect function of the currents (e.g. Eq.QI and is therefore 
also calculated at the center of the elements. In the next 
step, we need the curl of E, which naturally leads us back 
to the vertices and the results (dtB) immediately allow 
modifying B at those positions. Note that calculating 
the curl of E at the sample surface is not required. 

At the sample surface, the induction is obtained by 
integration, which is approximated by a summation over 
the discrete elements k within which J is constant. The 
three components (x, y, and z) of Eq. ([7]) may be written 
as 

H s , x (f 8 ) = £[j«4 k )(r s )- J z (k) 4 k) (r s )] (12) 

k 

H s ,y(r s ) = ^[J z (k) 4 k) (^)-4 k) 4 k) (^s)] (13) 

k 

H s Ars) = J2[ J i k) ^Hr s )-4 k) A^(f s )] (14) 

k 

The sums run over all elements of the grid, Jj denotes 
the i-component (i = x,y, or z) of the current density 

of the fc th element, and A; (r a ) the i-component of the 
function 



& k) ^) = T I rfV i 

4?r J (k) 



(15) 



where f s denotes a vertex position at the sample surface 
and the integral is performed over the volume of the k th 
element. We point out that it is essential to calculate 
these integrals numerically very accurately. Note, how- 
ever, that (r s ) depends only on the position vectors, 
hence the integrals need to be calculated just once dur- 



C. Magnetization loops of rectangular samples 

The model was successfully applied to fully 3D calcula- 
tions in conveniently short simulation times, as reported 
in Ref. 10. In the following, we focus on magnetization 
loops of rectangular samples in the presence of an homo- 
geneous external field, which occurs in most magnetome- 
try experiments (e.g. SQUID or VSM). The external field 
is applied in z direction, perpendicular to the top and 
bottom surface of our sample, and we assume isotropic 
material parameters (e.g. J c ) within the xy plane. In this 
case, the model may be simplified assuming the currents 
to flow exactly parallel to the nearest sample surfaces. 
The benefits of the procedure are still available in the 2 
dimensional calculations, making it possible to simulate 
realistic magnetization loops very rapidly (e.g. within 
about 1 to 10 minutes on a conventional standard PC). 
We point out that our assumption on the current flow 
is not rigorously valid (see Ref. Il8l ). but the deviations 
are only significant in very thin samples, which are not 
addressed here. 

For the actual implementation, we assume the sample 
center to coincide with the point of origin and the lat- 
eral surfaces to be parallel to the x or y axis. The fields 
and currents need to be calculated only at the cross sec- 
tion spanned by the z and y axis, on which a 2D grid 
is defined. Thus we apply almost the same procedure as 
illustrated in Fig. [3J but mapped on 2 dimensional el- 
ements (i.e. the cross section of each 3D element - see 
Fig. [3j. The calculation of the surface field (Eqs. [12]- 
[T4"]) is slightly modified to 



ff s ,x(r s ) - 

H a<y (f s ) = Y,J {m) Ci m) (r s ) 

m 

H s Ar s ) = $> (m) Ci m )(f s ) 



with 



^ (m) W) = — / d-v 

(m) 



,e*j x (r s - f) 



(16) 
(17) 

(18) 



(19) 



H s ^(f s ) vanishes due to the symmetry of the system. 
The sums run over all 2 dimensional elements (m) of the 
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FIG. 4: (Color online) Verification of the creep behavior in 
the simulations. Panel (a) present a magnetization loop (n 
= 20) followed by recording relaxation at 1 T for about 3 
hours. The inset shows the relaxation rate din m/d In i (~ 
— 1/n) together with a linear fit revealing n = 19.7. Panel 
(b) shows magnetization loops at different field sweep rates 
(10" 5 ,10" 4 ,10" 3 ,10 -2 , and 10 _1 T/s) and the inset the dy- 
namical relaxation rate din m/d ln/ioifa (— l/ n ) at IT, 
where the linear fit leads to n = 20.7. 



sample cross section. The integration is performed over 
one element (to) and - for the third dimension - over 
the closed current loop, exactly parallel to the nearest 
sample surfaces, ej indicates the unit vector of the a- 
priori known current flow direction. Thus C^ m \f s ) can 
also be accurately calculated during the initial stage and 
stored for later use. 



D. Tests of the simulations 

Several tests are available for verifying the output of 
the simulations. The first concerns H*, the lowest ap- 
plied field of a virgin m(H a ) curve, at which the cur- 
rents have penetrated the entire sample. It corresponds 
to the z-component of the field of that current distribu- 
tion at the center of the sample, which is easily calculated 
for constant currents (Bean model, e.g. Eq. [7j but with 
r = (0,0,0) T ). 

The simulations reveal H* by monitoring the current 
density at the sample center. The Bean model is approx- 
imated by applying constant currents and - to reduce re- 
laxation - a large n value. As expected, the deviation of 
H* from the direct calculation is reduced with increas- 
ing n value and almost vanishes (smaller than 1%) at 
n > 100. 

Secondly, we can test the creep behavior of the cur- 



rents by studying the relaxation of the magnetic moment. 
This time, the power law (Eq. 21) was applied with a 
field dependent critical current density and a constant 
n value (n = 20). Imitating the common experimental 
procedure, we started with simulating (parts of) a mag- 
netization loop, followed by recording the relaxation in 
the fully penetrated state, i.e. the time dependence of 
the magnetic moment - m(t) - at constant applied field 
(Fig. [4^) . According to theory the n value may be ex- 
tracted via | d In m/d In t | ~ l/n which matches very well 
(~ 1.5% deviation) the input value of n in the simulation 
(see inset of Fig. 2k) ■ 

A second way of studying relaxation is provided by the 
so-called dynamic relaxation, where magnetization loops 
are measured (or simulated) at different field sweep rates 
H a = dH a /dt (Fig. 2b)- The relaxation rate at constant 
applied field is now given by |dlnm/dln^o^a| — 
and is again found to be in good agreement (~ 3.5% 
deviation) with our input value of n (see inset of Fig. 2b). 

III. APPLICATION TO MAGNETOMETRY 
A. Evaluation of the critical current density 

1 . Experiment 

In the following we present two applications of the sim- 
ulations to common problems of magnetometry experi- 
ments. The first refers to the evaluation of the critical 
current density from magnetization loops. Usually the 
Bean model is applied, which presumably leads to certain 
deviations from the correct J c behavior. To analyze these 
deviations quantitatively we start from a given J C (B) be- 
havior and simulate magnetization loops. Then we apply 
the same evaluation procedure as used for the experimen- 
tal data to reevaluate J c , which is finally compared with 
our input J c curve. 

We start by explaining the experimental procedure 
used in most experiments for acquiring J c (77 a ) from mag- 
netometry measurements and introduce a simple exten- 
sion that allows to obtain J C (B). When determining the 
critical current density from SQUID or VSM measure- 
ments we usually do not know more than the magnetic 
moment as a function of applied field, m(_ff a ), and the 
sample dimensions a, 6, c, where we assume c to be paral- 
lel to H a . Thus, assumptions on J c are inevitable, which 
are in most cases a constant absolute value of J equal 
to J c and a flow direction parallel to the nearest sample 
surfaces (Bean model). The relation between m and J c 
is given by Eq. ([8]) leading to an analytic expression for 
rectangular samples in the fully penetrated stated that 
can be solved for J c 

jroil 4_ witna>6 (20) 

= abc and m\ denotes the irreversible magnetic mo- 
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ment generated by the critical currents, which is given 
by half of the hysteresis width of the magnetization loop. 
Equation (|20|) leads to a first approximation of J c as a 
function of the applied field - J c (H a ). To come closer 
to the actual material property, i.e., the critical current 
density as a function of the magnetic induction - J C (B) - 
(which is independent of sample geometry) , we addition- 
ally consider the field induced by the current distribution 
(cf. Ref. for the case of cylindrical samples) . Since the 
magnetic moment acquired by the magnetometry mea- 
surements always refers to the whole sample volume, the 
induction needs to be averaged in the same way to obtain 
m(B), i.e. 



with H s from Eq. ([7]), and 

(\H a + H s , z \) = - fd 3 r (fx e*j) z \H a + H s , z (f)\ (22) 
w J 

The integral runs over the whole sample volume, (fx e*j) z 
represents a weighting factor, which corresponds to that 
of J c in the calculation of m (cf. Eq. [8]) and reflects 
also the distance from the current element to the pick-up 
coils of the magnetometry devices, w is the integral of 
this factor over the sample volume. The weighting does 
not strongly affect the results, but is still significant for 
the behavior of J C (B). Note further, that we evaluate J c 
as a function of the absolute value of the B component 
parallel to the applied field (B z ). 

Finally, we address the evaluation of m\. The experi- 
ments provide m(H a ), which is usually the sum of mi and 
of the reversible magnetic moment m r . Additional mag- 
netic signals from the sample or the sample holder can 
usually be treated in a similar way as m T . mi is extracted 
from (half of) the difference of m at increasing (H+, B+) 
and decreasing (H-, £?_) fields at the same induction B. 
But note that the two branches of m(H a ) refer to differ- 
ent values of B at the same H a since the currents have 
opposite orientation in the two branches. Two situations 
arc distinguished. First, when the irreversible part domi- 
nates, as is the case in most experiments, B is separately 
calculated for each branch of m(H a ) by applying Eqs. ((20] 
- EH) witn 

mi ~ to. Then rrii(B) is obtained from 
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FIG. 5: (Color online) Results for different sample geometries: 
axixc=lxlx I mm 3 in (a), 3 x 3 x 3 mm 3 in (b), and 
1 x 0.5 x 0.1mm 3 in (c). The open symbols show the input 
J C (B) used to simulate the magnetization loops (e.g. inset of 
a), the dotted curve J C (H & ) and the solid line J C (B) evaluated 
from the magnetization curves by using Bean's model (see 
text). 



2. Simulation 



B = n (H a + H s J) (21) 



m(B-) - m(B+) 
mi(B) = (23) 

and the final J C (B) by evaluating Eq. (|20p again. In 
the second case, where the reversible parts are signif- 
icant or even dominate, the critical current density is 
usually low and therefore the corresponding field correc- 
tion (H SjZ ) small. Accordingly, we start with evaluating 
mi(H a ) — 0.5[m(iJ a _) — m(H a+ )], i.e. the irreversible 
magnetic moment from the hysteresis width at the same 
applied field, and then apply Eqs. (j2D1 - l2"Tj) for getting 
MB). 



The above experimental evaluation is not expected to 
result in the "true" J C (B) curve. The main reason for 
discrepancies is the (necessary) assumption of a constant 
current density within the sample, which is not fulfilled 
in realistic materials, since B may vary considerably and 
J C (B) is usually not constant. Accordingly, we expect 
more pronounced deviations, when the samples are large 
and have high J c , since both properties would enhance 
the (peak to peak) variation of B within the sample. At 
low fields this variation is comparable with (j,qH* , which 
may amount to several Tesla in typical "SQUID or VSM" 
samples. At higher fields, we expect better agreement 
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with the "true" J c , since J c becomes smaller and flatter 
as a function of B in most cases. 

To get some idea of the quantitative differences be- 
tween the "true" J c and J c from the above evaluation, 
we simulated a typical magnetometry measurement as- 
suming a sample size of 1 x 1 x 1 mm 3 , and a J c behavior 
described by 



Setting B = 0.2 T and a = 1/2 roughly leads to 
a behavior observed in many (melt textured) Y-123 
(YiBa2Cu307_2;) samples at low temperatures. J c q = 
1 x 10 10 Am~ 2 leads to [i H* ~ 2.1 T and a similarly 
large variation of B in the sample. 

Furthermore, we applied the power law (Eq. 3]) with 
E c = 5 x 10 _6 Vm _1 and n = 20 and a field sweep rate 
of /zo-ffa = 1 x 10 _2 Ts _1 . This rather fast field sweep 
rate corresponds to the maximum we can set in our VSM 
measurements and was chosen to reduce relaxation ef- 
fects, i.e. to get J ~ J c , where J is the real flowing 
current. 

The inset of Fig. [5^l presents the simulated magneti- 
zation curve of this sample up to fields of 9T, including 
the virgin curve. Before applying the methods of section 
IIII A 11 all parts of the curve, where the currents do not 
flow in the whole sample or do not have the same orien- 
tation need to be removed. This is done automatically 
by the evaluation program and includes the appropriate 
parts of the virgin curve (i.e. for H a < H*) and those 
parts directly after reversing the field, where the currents 
have opposite orientation. 

The results are illustrated in Fig. [5^. The open sym- 
bols show the original (input) J C (B), i.e. the "true" 
J C (B), according to Eq. The dotted line presents 

J°(H a ), where the subscript "e" indicates that J c has 
been evaluated from the magnetization curve (of the in- 
set of Fig. [SJt) with the (experimental) methods of section 
IIII A 11 and the solid line is the final result Jc(B) when 
considering also the field correction according to Eq. (|2"Tj) . 
Jc(H a ) shows significant deviations from J C (B) at fields 
H a < H*. J£(B), however, matches quite well despite 
the rather rough assumption of a constant current den- 
sity in the sample (note that the peak to peak variation 
of B in the sample reaches 2.4 T at constant applied field 
in this specific simulation). Slower field sweep rates or 
field-step measurements would shift the curves to slightly 
smaller values due to relaxation. Note further that the 
evaluation method does not give access to J c at or close 
to B = 0, since the mean induction (cf. Eq.l2"2"]) is always 
larger than zero. 

To get further insight, we repeated the procedure for 
different samples and J C {B) curves (e.g. with exponential 
behavior) and found similarly good agreement as above. 
For instance, Fig. [5]b shows results on a sample with di- 
mensions 3 x 3 x 3mm 3 and the same J C (B) as in the 
above example (only E c is slightly adapted to account 
for the larger sample size). Although hqH* of 4.4 T is 



quite large and J^{H a ) starts to deviate from J C (B) at 
a similarly high field, these deviations are again nicely 
corrected by Jc(B). 

Analyzing situations with much larger fJ-oH* does not 
make sense, since the evaluation of J c requires external 
fields larger than /iqH*, but typical experimental devices 
do not provide fields much above 5 - 10 T. Increasing 
J C (B) instead of the sample size leads to similar effects. 

As a final example, we show results on a sample with 
a size of 1 x 0.5 x 0.1mm 3 , which is representative for 
a typical single crystal geometry. We took the same 
J C (B) behavior as in the above examples but with J c o 
= 5 x 10 10 Am- 2 , accordingly fi H* ~ 2.1 T. Fig. [Sfc 
presents the results, which again demonstrate the ex- 
cellent agreement between the evaluated J°(B) and the 
input curve J C (B). 



B. Superconducting volume fraction 

The evaluation of the critical current density as de- 
scribed in the previous sections is only valid if the cur- 
rent flow is unimpeded over the whole sample. This 
is illustrated by Eq. ([20l . which shows that a sample 
with dimensions a x b x c would induce a larger mag- 
netic moment than the sum of two samples with dimen- 
sions a x 6/2 x c, with the same J c as the large sample. 
Thus, grain boundaries and other macroscopic inhomo- 
geneities (e.g. normal conducting inclusions), which im- 
pede the current flow, would lead to an underestimation 
of J c , when applying the above method. A similar effect 
is caused by overestimating the superconducting sample 
volume, e.g. when the material properties are degraded 
at the surface. We further note, that magnetization mea- 
surements are not only used for evaluating J c , but also 
for determining all kinds of superconducting parameters, 
such as the critical magnetic fields or the characteristic 
lengths, etc. (e.g. Ref. [2l|). Most of these quanti- 
ties are derived from the magnetization, M = m/Q, i.e., 
it is again important to exactly know the volume (f2). 
We point out that both surface degradation and grain 
boundaries were even found in samples believed to be 
single crystals. In the following we show that the reverse 
branch of a magnetization loop is very sensitive to such 
imperfections, and that we can prove whether or not a 
sample is single grained, and even detect, if parts of the 
surface are degraded. 

We started our analysis with experimental magnetiza- 
tion loops of a Y-123 sample at 30 K. The size of the rect- 
angular sample was precisely measured with an optical 
microscope and found to be a ~ 1.44 mm, b ~ 0.88 mm, 
and c ~ 0.76mm. Using these data, J C {B) was calcu- 
lated with the above method (section IlII A lj) . and the 
result was used for simulating the magnetization loops. 
As expected, there is agreement with experiment over 
most part of the field range. Small differences are only 
observed near H a = 0, since J c close to B = needs 
to be extrapolated. We are, however, not interested in 
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FIG. 6: (Color online) Comparison of the magnetic moment 
from experiment and simulations with different sample surface 
sizes. The arrows indicate the interval, within which currents 
with opposite orientation flow in the sample. 



this region, but look closer to that part of the reverse 
branch of the loop, where currents with opposite orien- 
tation flow in the sample (indicated by the two arrows 
in Fig. [6]) and which had been disregarded for evaluat- 
ing J c . As shown in Fig. [5] there is perfect agreement 
between experiment (open circles) and simulation (solid 
line, a x b = 1.44 x 0.88 mm) even in this part of the 
loop, thus demonstrating that the size of the sample sur- 
face had been determined correctly and that the sample 
is pretty homogeneous for the macroscopic current flow. 
Comparing only the magnetization slope directly after 
the reversal point (up-arrow in Fig. [6]) (inset of Fig. [6]) 
provides information on the sample size and surface, but 
not on the interior since the newly penetrating currents 
(which have opposite orientation with respect to the ex- 
isting currents) flow only at the sample surface. With de- 
creasing applied held, the new currents penetrate into the 
sample and as soon as the current front reaches a region 
with a different J c (e.g. at a grain boundary or a nor- 
mal conducting inclusion) the macroscopic current path 
and therefore also the magnetic moment would start to 
deviate from the simulations (where perfect homogeneity 
is assumed). In addition, the field, where the new cur- 
rents arrive at the center (down-arrow in Fig. O, would 
be different. To conclude, only comparing experiment 
and simulation over the whole penetration interval (i.e. 
form the up- to the down-arrow in Fig. [6]) can clarify if 
samples are homogeneous (at least for the current flow) 
over the whole volume. 

Let us assume that the measurement of the sam- 
ple size had resulted in slightly larger values, namely 



a = 1.49 mm and b = 0.93 mm, i.e. only 0.05 mm for 
both lengths, which increases the volume by about 9%. 
We repeated the whole procedure with these data (which 
leads to a slightly smaller J c ) and find again agreement 
on most parts of the m(H a ) loops. However, inspecting 
the reverse branch more closely, we find significant dis- 
agreement between the two arrows (Fig. [5]), i.e. where 
not all currents are equally oriented, which demonstrates 
different sample sizes. This indicates either an error in 
the measurement of the sample size, or the fact that small 
regions of the lateral surfaces are non-superconducting. 
While a smaller J c hardly affects the slope of m(H a ) close 
to the reversal point, the larger sample size and the cor- 
respondingly larger current loop enhances the (absolute 
value of the) slope by about 10% as shown in the inset 
of Fig. [HI (note that H* also increases). 

Secondly, we studied the case of a sample divided into 
two independent grains (bicrystal) with sizes 0.72 x 0.88 x 
0.76mm 3 (i.e. a/2 x b x c). Since the relative contribu- 
tion of the critical currents to the magnetic field is small 
at high applied field, any interaction between the two 
grains (i.e. magnetic fields from the other grain) was 
neglected. Having the same overall magnetic moment 
(which is simply the sum of m from the two grains) leads 
to a larger J c than in the large (single grained) sample. 
As expected, the slope of m(H a ) at the reversal point is 
significantly flatter in the bicrystal (about 17%) than in 
the single grain sample (inset of Fig. [fj]) making it easy 
to distinguish between the different situations. 

Finally, we compare our results on the slope of the 
reverse branch with the well known analytical expres- 
sion given by Angadi et al22, i.e. dm/diJ = — tt 2 R 3 /Q, 
with = In (8R/c) — 0.5 and R the radius of a cylinder 
which has the same surface area as the cuboid. Using 
this method results in a slope that deviates from simu- 
lation and experiment by only about 5% (inset of Fig. 
|B]). However, applying the same method to the bicrystal 
leads to a slope that actually matches the exact result of 
the single grain sample even better. Thus, this approxi- 
mation cannot be recommended to verify if samples are 
single grained or decomposed into a few grains. 



IV. SUMMARY 

Details of a simple method for two or three dimensional 
simulations of the current dynamics in a superconduc- 
tor have been presented. Solving the differential form 
of Maxwell's equations in the sample interior makes the 
calculations very fast. The integral form is only used at 
the sample surfaces to meet the correct boundary con- 
ditions. The procedure results in the time and position 
dependence of the magnetic induction and the critical 
current density, which makes it very convenient to ap- 
ply field dependent material parameters, like J C {B). The 
method allows simulating different kinds of experiments, 
like magnetization loops, static and dynamic relaxation 
measurements, etc., in which also the reversible super- 
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conducting properties may be considered. 

To demonstrate the potential of the method, we pre- 
sented two examples, which are important for standard 
magnetometry measurements. In the first, magnetiza- 
tion loops were simulated by applying a given J C (B) and 
then evaluated by the methods usually applied to exper- 
imental data, i.e. the Bean model. Comparison of the 
given J C (B) and the 'Bean' J C (B) shows good agreement 
even for large samples and high J c , if the 'Bean' result is 
corrected by a simple numerical expression. 

Secondly, we showed that surface degradations, which 
reduce the superconducting volume fraction of a sample, 
or inhomogeneities like grain boundaries, which influence 



the current flow, can be detected by carefully comparing 
the magnetization loops from experiment and simulation 
at regions, where currents with opposite orientation flow 
in the sample. 
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